DYNAMIC BEHAVIOR OF STOCHASTIC GENE EXPRESSION 
MODELS IN THE PRESENCE OF BURSTING 



M. C. MACKEYt, M. TYRAN-KAMINSKA* , AND R. YVINEC§ 



Abstract. This paper considers the behavior of discrete and continuous mathematical models 
for gene expression in the presence of transcriptional/translational bursting. We treat this problem 
in generality with respect to the distribution of the burst size as well as the frequency of bursting, 
and our results are applicable to both inducible and rcprcssiblc expression patterns in prokaryotcs 
and eukaryotes. We have given numerous examples of the applicability of our results, especially in 
the experimentally observed situation that burst size is geometrically or exponentially distributed. 
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1. Introduction. Recent spectacular advances in the ability of experimentalists 
to monitor the temporal behavior of single molecules [U[E1[18 ( 22, 23,|2H1|33] inside cells 
has led to a quantum leap in our knowledge of their behavior as well as a plethora 
of data that challenge mathematicians. These techniques are so refined that they 
allow the single molecule quantification of the transcription of mRNA as well as the 
translation of the mRNA into protein. This visualization has shown that in many 
cases these transcription and translation processes occur in quantal bursts in which 
a few molecules are produced during a discrete period of time. An analysis of the 
data obtained from such experiments has given us many details of the nature of 
the bursting kinetics that are being used to guide mathematical modeling of these 
fascinating processes. 

This paper utilizes the two main approaches that have been employed to model 
these bursting processes, i.e. a discrete formulation for the numbers of molecules [27] 
or a continuous one [6] I15j and illustrates the common features of both as well as the 
differences. Modeling (as opposed to simulation [16] which we do not consider) of 
the details of gene expression as a discrete Markov process has an extensive literature 
(c.f P HUJ 1201 HH HS1 12U EH ED]) that has recently seen a flurry of activity. The 
other approach that has received extensive attention is modeling of the process as a 
continuous one and [5] [T3] [121 El 12S] are representative of these efforts. The reader 
can consult for an excellent expository account of the connection between these 
two approaches. 

In the discrete Markov models, steady-state analytical solutions of the master 
equation can often be obtained using the moment generating function. For the contin- 
uous model formulations, one needs to solve the Fokker-Planck-like equations, some- 
times using Laplace transforms. When solutions are not available, moment equations 
can be derived and usually solved [TH [29]. Though continuous models have many 
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analytic advantages over discrete ones, it is also the case that information of poten- 
tial importance may be lost in the continuous model formulation compared with the 
discrete formulation. 

This paper presents a general one dimensional model for bursting gene expression 
in both a discrete Markov process formulation as well as a continuous situation. Sec- 
tion[2]prcsents some general background material while Section |3] presents the discrete 
version of the bursting model. Section [3 . 1 1 develops the general formulation of the dis- 
crete model while Scction l3~2"l deals with the special case in which the burst amplitudes 
are geometrically distributed. Section 0] develops the corresponding continuous model 
of the bursting expression, with a general development in Section |4~T1 and Section l4~2l 
devoted to the situation where the burst amplitudes are exponentially distributed-a 
situation often found experimentally. Section 14.31 concludes with an examination of 
a generalization of the exponential distribution of burst amplitudes. The paper ends 
with some general observations in Section [5] Throughout the paper, our results are 
illustrated with numerous examples. 

2. Notation and background. Let the triple (E,£,m) be a cr-finite measure 
space and let L 1 = L 1 (£',f , m) with norm denoted by || ■ ||i. A linear operator P on 
L 1 is called sub stochastic (stochastic) if Pu > and ||Pu||i < ||u||i (||Pu||i = ||u||i) 
for all u > 0, u G L . We denote by D the set of all probability densities on E, i.e. 

D = {u G L 1 : u > 0, = 1}, 

so that a stochastic operator transforms a density into a density. In the particular case 
of a countable set E with £ being the family of all subsets of E and m the counting 
measure, the space L 1 will be denoted by t 1 . 

Let V : E x £ — > [0, 1] be a stochastic transition kernel, i.e. P(x, •) is a probability 
measure for each x G E and the function x i— > V(x, B) is measurable for each B G £, 
and let P be a stochastic operator on L 1 . If 

/ Pu(x)m(dx) = / V{y,B)u(y)m(dy) for all B e £, u S D, 
J b J E 

then P is called the transition operator corresponding to V . A stochastic operator 
P on L 1 is called partially integral or partially kernel if there exists a measurable 
function p: E x E — > [0, oo) such that 

p(x,y) m(dy) m(dx) > and Pu(x) > / p(x,y)u(y) m(dy) 
'eJe Je 

for every density u. If, additionally, 

p(x,y)m(dx) = 1, y G E, 

E 

then P corresponds to the stochastic kernel 

V(y 7 B)= [ p(x 1 y)m{dx) 1 V eE,Be£, 

JB 

and we simply say that P has kernel p. Note that each stochastic operator on t 1 has 
a kernel. 

We denote by T> (A) the domain of a linear operator A. We say that A C B, 
or that B is an extension of A, if V(A) C V(B) and Bu = Au for u G V(A). The 
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operator A is said to be closable if it has a closed extension. If A is closable, then the 
closure A of A is the minimal closed extension of A; more specifically, it is the closed 
operator whose graph is the closure in L 1 x L 1 of the graph of A. For an exposition 
of semigroup theory we refer to [5]. 

A semigroup {P(t)} t > of linear operators on L 1 is called sub stochastic {stochas- 
tic) if it is strongly continuous and for each t > the operator P(t) is substochastic 
(stochastic). A density u* is called invariant or stationary for {P(t)}t>o if u* is a 
fixed point of each operator P(t), P(t)u* = u* for every t > 0. 

Theorem 1 ( [211 Theorem 2]). Let {P(t)} t >o be a stochastic semigroup such that 
for some to > the operator P(fo) is partially integral. If the semigroup {P(t)} t >o 
has only one invariant density u* and u* > a.e. then 

lim \\P{t)u - u* ||i = for all u£D. 



3. A discrete bursting model formulated as a Markov process. This 
section considers bursting gene expression as a Markov process. 

3.1. The general case. In this section we model the number of gene products as 
a pure-jump Markov process X = {X(t)} t >o in the state space E = {0, 1,2,.. .}. Thus 
a master equation governs the dynamics evolution of probabilities. A general one- 
dimensional bursting gene expression model |llj may be constructed as follows. Let 
n be the number of gene products and P n (t) = Pr(X(t) = n) denote the probability 
of finding n gene products inside the cell at a given time t. We shall include a loss 
(n — > n — 1) and gain (n — > n + k) of functional processes in terms of the general rates 
7„ and A n , respectively. The step size assumes the values k = 1, 2, . . . and is a random 
variable (independent of the actual number of gene products) with probability density 
function h, so that X/fc^i hk = 1- Therefore, our general master equation describing 
the time evolution of the probabilities P n to have n gene products in a cell is an 
infinite set of differential equations 

dP " 

(3.1) = Jn+lPn+l - InPn + ^ /lfcA n _feP„_fc - A„P„, 71 = 0,1,..., 

fc=l 

where we use the convention that Y2k=i = ®- ^ c supplement f|3 . 1 [) with the initial 
condition P„(0) = v n , n = 0, 1, . . ., where v = (i>„) n >o € I 1 is a probability density 
function of the initial amount A(0) of the gene product. In the following paragraphs, 
we consider the existence and uniqueness of solutions of (|3.ip together with conver- 
gence to a stationary distribution and then summarize our results in Theorem [2j 
Assume that 

+ oo 

(3.2) A o >0, 70 = 0, 7«>0, A„,/i n >0, n = l,2,..., ^/i n = l. 

n=l 

The process X is the minimal pure jump Markov process with the jump rate function 
tp(n) = \ n + 7„, n > 0, and the jump transition kernel /C given by 

{q n , if j = -l,n> 1, 

{l-q n )hj, ifj>l,n>0, q n = . 
0, otherwise, An + 7,1 
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First, we recall the construction of X. Let {£, k } k >o, be a discrete time Markov chain 
in the state space E = Z + = {0, 1, . . .} with transition kernel /C and let {e k }k>i be 
a sequence of independent random variables, exponentially distributed with mean 1. 
Set To — and define recursively the times of jumps of X as 

T fc = T fc _ 1 + £k fc = l,2,.... 

<A?fc-i) 

Starting from X(0) = £o we have 

X(t) = £ k , T k <t<T k+1 , k = 0,1,2,..., 
so that the process is uniquely determined for all t < Too, where 

^ = lim T k 

k^roo 

is called the explosion time. If the explosion time is finite, we can add the point — 1 
to the state space and we can set X(t) = —1 for t > Too- The process X is called 
nonexplosive if Pi (Too = oo) = 1 for all i £ E, where Pi is the law of the process 
starting from X(0) — i. In particular, if the chain {£fc}A;>o is recurrent, then X is 
nonexplosive. 

We now rewrite equation (|3.1|) as an abstract Cauchy problem in the space I 1 . We 
make use of the results from [32]. Let K be the transition operator on 1 1 corresponding 
to K defined as in (|3.3|) . For v = (v n ) n >o £ £ l we have (Kv)q = q±vi and 

(Kv) n = q n +iv n+ i + ^2 h k(l - q n -k)v n -k, n = 1, 2, 

fc=l 

Define the operator 

oo 

Gu = —ipu + K(ipu) for u £ l 1 ^ = {u £ I 1 : y n |wn| < oo}. 

n=0 

There is a substochastic semigroup {P(t)} t >o on I 1 such that for each initial proba- 
bility density function v £ the equation 

nil 

(3.4) ^ =G(u) ' <>0 ' u (°) = u ' 

has a nonnegative solution u(t) which is given by u(t) = P(t)v for t > and 

oo 

(P(t)v) n =^P 3 -(l(t) = n,t < Tojvj, n = 0, 1, ... . 

3=0 

The process X is nonexplosive if and only if the semigroup {P(t)}t>o is stochastic. 
Equivalently, the generator of the semigroup {P(t)} t >o is the closure of (G,£^). In 
that case the solution u(t) of (|3.4[) is unique and it is a probability density function 
for each t, if v is has these properties. 

The equation for the steady state p* = {Pn)n>o of (|3 . 1 [) is of the form 

n 

(3.5) Jn+lP^+i ~ lnP* n + ^ h k^n-kP* n - k ~ KPn = 0, 71 = 0,1,.... 

fc=l 
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Observe that Jipl = AoPg an( i that we can rewrite (|3.5[) as 

n-1 

7n+U>n+l _ 7nK = A «-Pr* - X! h n-k^kP* k , H = 1, 2. 

Hence 

n n j — 1 

7-i+iP,*+i = XI A J'Pi ~ X E h J-k X kPt 

3=0 3=1 fc=0 

and changing the order of summation, we obtain 

1 " - 

(3.6) p* +1 = ^— — ^Thn-k^kPk, n = 0,l,.. 

where 



foi = X /ij, Z > 0. 

j=l+l 

Thus given Pq, equation p.6j) uniquely determines p*. Consequently, there is one, and 
up to a multiplicative constant only one, solution of equation (|3.5[) . If j»q > and 
either hi > for all / > or A; > for all Z > 1, then p* n > for all n > 1. Now, if 

CJO oo 

(3.7) ^Pn = 1 and X]( A " +7«)K < °°> 

n=0 n=0 

then p* G G(p*) = 0, and K(ipp*) = i/jp*, which implies that the semigroup 
{P(t)} t >o is stochastic. We have thus proved the following result, which is an analog 
of Theorem [1] for the discrete bursting model. 

Theorem 2. Assume condition (|3.2[) and suppose that a strictly positive p* = 
(Pn)n>o given by (|3.6p satisfies (|3.7[) . Then for each initial probability density function 
v = (v n )n>0 S equation (|3.ip /ias a unique solution which is a probability density 
function for each t > and satisfies 

oo 

£— >oo z — ^ 

71 = 

Remark 1. From (13.61) it follows that 



E^i = EE E fciUrf = EE E 

n=0 n=0fe=0 \7=n-fc+l / k=0n=k \j=n— k+1 

The mean value E(/i) of the distribution h can be represented as 

oo oo oo 

E(/o = £>i = E E v 

j— n— j— n+1 
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Wc thus obtain 



oo oo 



E w = E E h i 

n—k j—n — k-\-l 

for each k > 0. Combining these leads to 

oo oo 

E 7n+lK+l = E ( h ) E XkP *k- 
n=0 k=0 

3.2. Bursting with a geometric distribution. Next, we give sufficient con- 
ditions for (|3.7p in the case when h is geometric 

(3.8) h k = (1 - b)b k -\ k = l,2,..., 

with b e (0, 1). Since 



E ^ = * 

j— n— k+1 

we obtain the following equation for p* = (Pn)n> 



-k 



(3.9) — — = , n = 0,l.... 

Pn 7n+l 

Explicit stationary solutions in this case were recently obtained in pQ. However, for h 
geometric, we can go further and prove convergence to this stationary state with the 
following result which follows from Theorem [2] and Remark [TJ 

Corollary 3. Assume that condition (|3.2j) holds. Suppose that h is geometric 
as in (|3.8|) . Then p* = (p*)„>o is given by 

,o 1n x * * TT A fc-i + b ^-i Po x o 'fx 1 A fc + bj k 

(3.10) P„=PoM = II , n=l,2,.... 

In particular, if 

(3.11) lira sup — < 1 — and liminf7„ > 0, 

£/ien £/ie conclusions of Theorem [H ZioZd. 

Example 1. Consider X n to be a Hill function of the form 

1 + Bn N 

< 3 ' 12 > A " " 

where A, A, N > and > 0. If h is geometric and 

A9 



lim inf 7„ > 



A(l-6)' 



then condition p. lip holds. 

Remark 2 (^Bifurcation in the discrete case ). Equation l3.9l can be used to examine 
the bifurcations in the stationary density, defined as changes in the number of maxima, 
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as a function of the model parameters. The number of maxima are linked to the 
number of sign changes of 

(3.13) n h-> A„ + b~f n - 7„ + i. 

In particular, p* has a maximum at if Ao < 7 i, and each successive sign change of 
(|3.13p gives a maximum/minimum of p* . 

We now provide examples for which the stationary distribution can be identified 
explicitly. In the following examples we assume that h is geometric with parameter b 
as in (|3.8[) and that j n = jn, n > 1, with 7 > 0. 

Example 2 (Negative binomial ). Suppose that A„ = Ao + An with Ao > 0, A > 0. 
We have A„ > for each n. Substituting 7^ and Afe into (|3.10p gives 

* *Stt7 Ao 4.^ f x+b A n n 1 

^ = ^{^Tx + k ){—) ' n = ' 1 '- 

Thus p* G t 1 if and only if 

A + fry < 7. 

In that case we obtain the negative binomial distribution 



K=(^V(i-p) a , n = ,i, 



n\ 
where 

A + 67 A 



V 



7 ' 67 + A ' 

and (a)n is the Pochhammer symbol defined by 
Via + n) 

(a)n = w x = o(a + l)(a + 2) . . . (o + n - 1), (a) = 1. 

r(o) 

This was previously obtained in |27j . 

Example 3 (Mixture of logarithmic distribution ). Suppose that Ao > and A„ = 
for n > 1. Then 

* * A b"^ 1 
Pn=Po ) "- = 1,2,..., 

7 n 

which can be rewritten as 

b n 67 
P » = - nWl-b) {1 - Po) ' n=1 > 2 -- P ^ 6 7 -A ln(l-6) - 

The distribution 

b n 

pa = 0, p n = — — , n = l,2,..., 

n ln(l — 0) 

is called a logarithmic distribution. 
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If we assume that A„ = for n > m, then we obtain the following distribution 
Pn=Po^Ii(j^ + k ) > n = 0,...,m, 



and 



rn \ * — * 



J 

j=0 



n > m, 



where c and Pq are such that 



j=m+l ■' j=0 



In particular, this type of distribution will be obtained if we take Ao > 0, A < 0, and 

An. = 



Ao + An, if n < — Ao/A, 
0, otherwise. 

Example 4 (Hypergeometric distributions/ We now take 

where A > 0, A > 1, 9 > A. We find that, for each n, 

A 7l + 6771 b(n + a 1 )(n + a 2 ) 



where 



and 



7 n + 6i 



hi = ^, ai = i (a - /3) , a 2 = ^(a + /3) 



A A9 2 4A 

/r = a - 



A 67A 67 A 

Since A > 1 and 9 > A, we can find a nonnegative j3, thus a 2 > ai > 0. Consequently, 
the stationary distribution is of the form 

3.15 p = _— — — - — — n = 0,l,..., 

2 Fi(ai,a 2 ;6i;6) (&i)„ n! 

where 2 i*i is Gauss' hypergeometric function 

{a\) n (a2)n x n 



2 Fi(ai, a 2 ; 61; a;) = 



(&i)„ n! 



Example 5 (Generalized hypergeometric distributions/ The generalized hyper- 
geometric function p F q is defined to be the real analytic function on R given by the 
series expansion 

p F q ( ai , . . . , a p ; bi, . . . , b q , x) - ^ {hi)n ^ -j-. 
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The negative binomial distribution in Example [5] for the case of A = has the probabil- 
ity generating function s h-> iF (ai;bs)/ iF (ai;b) with ai — A0/&7. The distribution 
obtained in Example 0] has the probability generating function 

2Fi(ai,a 2 ;bi;bs) 
2 Fi(ai,a 2 ;bi;b) ' 

Extending both of these examples we suppose that A„ > is a rational function of n 
satisfying 

A» + bjn _ (re + ai) . . . (n + a g+ i)6 _ 
7 (n + 6i)...(n + 6 g ) 

Then p* = (p* )„>o has the probability generating function 

9+ iF 9 (a 1; . . . , a 9+ i; 61, . . . , 6 g ; 6s) 



9+1^9 



F q (a!, . . . , a q+ i;bi, . . . ,b q ;b) 



4. Continuous bursting model. 

4.1. The general case. In this section we consider a continuous state space 
version of the model presented in Section[3j which is a piecewise deterministic Markov 
process (PDMP) Y = {Y(t)} t >o with values in E = (0, 00) where Y(t) denotes the 
amount of the gene product in a cell at time t,t > 0. We assume that protein molecules 
undergo degradation at a rate 7 that is interrupted by production at random times 

h<t 2 < ... 

occurring with intensity ip, and that both ip and 7 depend on the current number of 
molecules. At each i& a random amount of protein molecules is produced, so that the 
process changes from Y(tk— ) to Y(tk) = Y(tk~) + efe, k = 1, 2, . . ., where {efe}fe>i is 
a sequence of random variables such that 

Pr(e fc G B\Y(t k -) = y) = [ h(x, y)dx, 

JB 

where ft. is a nonnegative measurable function satisfying 

POO 

(4.1) / h{x,y)dx = l, y>0. 

Jo 

The time-dependent probability density function u(t, x) is described by the continuous 
analog of the master equation P31 [T5] 

(4 2) du(±x± = d(j(x)u(t,x)) _ ^ {x)u ^ x) + ^ y)h{x _ ^ y)dy 

with the initial probability density u(0,x) = v(x), x > 0. 
We assume that 7 is a continuous function such that 

r s dx 

(4.3) 7(x) > for x > 0, / — - = +00, 

Jo 7W 
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for some S > and that ip is a nonncgative measurable function with ip/j being 
locally integrable on (0, oo) and satisfying 

(4.4) I" ^\dx = +oo. 

Jo 70) 

From (|4.3|) it follows that the differential equation 

(4.5) x'(t) = -7(z(f)), i(0) = a: > 0, 

has a unique solution which we denote by7r t a;,i>0,3;>0. For each x > we have 
ir t x > for all t > and 7r t x — > as i — ► oc. This and condition (|4.4|) give 

/ (pOjr s x)ds = / — j — -dy — > oo, as t — > oo, 
Jo A t x 7(2/) 

which implies that the function 

t I — y 1 — e -/o ^(Ta^ds 

is a distribution function of a positive and finite random variable for every x > 0. 

We now recall the construction of the minimal piccewisc deterministic Markov 
process Y (see e.g. [3J[3] or [35] for details). Let {£fc}fc>i be a sequence of independent 
random variables exponentially distributed with mean 1, which is also independent of 
{efc}fc>i- Set to = 0. For each k = 1, 2, . . . and given Y(t k -i) the process evolves as 

(A(!\ V f f \ _ / 7rf-tfc_l^(*A!-l), t k -i<t<t k , 

( ' \ Y(t k -) + e k , t = t k , 

where t k = t k -i + At k and At k is a random variable such that 

Pr(Ai fe < t\Y(t k -i) = x) = l-e-fX"'*')*', t,x>0. 

The random variable At k can be defined with the help of the exponentially distributed 
random variable e k through the equality in distribution 

£ fe = / tp(-K s Y{t k -i))ds, 



JO 

which can be rewritten as 

e k = Q(n Atk Y{t k ^)) - Q(Y(i fc _i)), 
where the non-increasing function Q is given by 

(4.7) Q(x) = [ ^±dy, 

and x = +00, when the integral is finite or any x > otherwise. Since Y(t k —) = 
^At k Y(tk—i), we obtain the following stochastic recurrence equation for {Y(t k )} k >o 



(4.8) Y(t k ) = Q- 1 (Q{Y{t k . 1 ))+s k ) + e k , fc = l,2,..., 
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where Q _1 is the generalized inverse of Q, Q 1 (r) = sup{x : Q{x) > r}. Consequently, 
Y(t) is defined by (|4.6[) for all t < t^, where = lim/^oo t^ is the explosion time. 
As in the discrete state space we can extend the state space E by adding the point 
— 1 and define Y(t) — — 1 for t > too. Let be the law of the process Y starting at 
Y(0) = x and denote by K x the expectation with respect to P x . 

Remark 3. Note that if condition (|4.4j) holds (cquivalently Q(0) = oo) then the 
amount of the gene product {Y(tk)}k>o at the jump times is a discrete time Markov 
process with transition probability function given by 



K(y, B)= [ k(x, y)dx, B e B((0, oo)), 

JB 



where 



(4.9) k(x,y)=eQM [ V l (0) ( z )h(x-z,z)^-e-^dz, x,y>0. 

Jo i{ z ) 

If Q(0) < oo then the random variable Aii is infinite with positive probability, since 
we have for any x > 

Pr(A*i = oo|Y(0) = a;) = lim Pr(A<i > t\Y(0) = x) = e Q(*)-Q(°) > o, 

t— >oo 

which then forces the process Y(t,u) starting form Y(Q,oo) = x to be 7r t (cc) for all t, 
if lu is such that Ati(o>) = oo. 

In what follows we assume that (|4.3[) and (|4.4|) hold. We rewrite equation (|4.2p 
as an abstract Cauchy problem in L 1 

nn 

(4.10) Tt =Cu > u ^ = v > 

where the operator 

(4.11) CU{X) = ^^{X)) _ v{x)u{x) + j* v{y)u{ ^ y)h{x _ y } y)d y 

is defined on the domain 

(4.12) V(C) = {u e L 1 : 7 u e AC, ( 7 u)' G L 1 , lim (j(x)u(x)) = 0, ipu £ L 1 }, 

and ju S AC means that the function x i— > 7(x)u(a;) is absolutely continuous. 
From [T3J [32] it follows that there is a substochastic semigroup {P(£)} t >o on L 1 
such that for each initial density v £ T>(C) equation (|4.10p has a nonnegative solution 
u(t) which is given by u(t) = P(t)v for t > and 

(4.13) f ¥ x (Y(t) e B,t <t OQ )v(x)dx = [ P(t)v(x)dx 
Jo Jb 

for all Borcl subsets B of (0, oo). The semigroup {P(t)} t >o is stochastic if and only 
if its generator (C,D(C)) is the closure of the operator (C, D{C)). 

We first study the fixed points of the semigroup, showing that {P(t)}t>o has no 
more that one invariant density through 

Proposition 4. The substochastic semigroup {P{t)}t>o can have at most one 
invariant density. 
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Proof. Recall that u* is an invariant density for the semigroup {P(t)} t > if and 
only if it is an invariant density for the resolvent operator 

poo 

Rv := R(l, C)v = / e-tp^vdt. 
Jo 

The operator R is substochastic and it satisfies Rv > R\v for any nonncgative v 6 L 1 
(see [H]), where 

i r°° 

R lV (x) = -4t / v{y)e Q{y) - Q{x)+I * ^ dz dy, x > 0. 

Note that i?i is the resolvent operator R(l, A) of a substochastic semigroup {S^)}^ 
with generator 

. , . dH(x)u(x)) . . . , „,„n 
4u i = v , v ;; - <p(x)u(x), u e D C . 
dx 

Since for any two nonncgative and nonzero Vx, u 2 £ we can find c(v.i) > such that 

Vi(y)dy>0, i = 1,2, 

c(«0 

we obtain Rvi(x) > for all x < min{c(vi), 0(^2)}, i = 1,2. Now suppose that ui,U2 
are densities such that u = u\ — u 2 is nonzero. Then both u + = max{0, u} and u~ = 
max{0, — u} are nonnegative and nonzero. Thus, R(u + )(x) > and R(u~)(x) > for 
x < c and some c > 0. We have 

\Ru(x)\ = \R(u + )(x) - R(u~)(x)\ < R(u+)(x) + R( U -)(x) = R(\u\)(x), 

thus the inequality is strict on a set of positive measure, which implies that if u\—u 2 7^ 
then 

\\Rui - Ru 2 \\i < \\R\ui - u 2 \\\i < \\ui - ua||i- 

Consequently, the operator R can have at most one invariant density. □ 
Let K be the transition operator on L 1 given by 

/•OO 

(4.14) Kv(x) = I k(x,y)v(y)dy, v G L 1 , 

Jo 

where the kernel k is as in (|4.9[) . Observe that 

(4.15) Kv(x) = [ h(x- 2,z)44e" Q(z) f°° v{y)e Q{y) dydz. 

JO l( z ) Jz 

A mild condition on the transition operator K, in conjunction with Theorems 3.6 and 
5.2 of [35], has interesting consequences for {P(t)}t>o as contained in the following 
result. 

Proposition 5. If the transition operator K is mean ergodic, i.e. for any v G L , 
v > the sequence 



1 n— 1 

- V K*v 
n ^ 



3=0 
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is convergent in L 1 , then the semigroup {P(i)} t > is stochastic. 

In particular, if K has a strictly positive fixed point, i.e. there is v* such that 
Kv* = v* and v* > a.e., then K is mean ergodic [12]. Note that a mean ergodic 
stochastic operator has a nonzero fixed point. 

We now describe invariant densities for the semigroup {P(t)} t >o with the help of 
fixed points of the operator K . 

Theorem 6. Let 

/•OO 

H{x,y) = / h(z,y)dz, x > 0. 

J X 

Suppose that there is a nonnegative solution u* of the equation 

(4.16) j(x)u*(x) = H(x - y,y)ip(y)u*(y)dy 

Ja 

such that ipu* G L . Then the function 

(4.17) v*(x)= [ h(x ~y,y)(f(y)u*(y)dy 



is a fixed point of the operator K in L , where K is as in (|4.15p . Moreover, ifu* G L 1 
then u* G T)(C) and C{u*) = 0, where C is as in (|4.1ip . 

Conversely, if the operator K has a nonnegative fixed point v* G L 1 then the 
function 

(4.18) u*(x) := — - / e Q(y) - Q(x) v*(y)dy 

7W Jx 

is a solution of (|4. 16[) and ipu* G L 1 . 

Proof. Let u* be a solution of (|4.16[) such that ipu* G L 1 . Since 

lim h y oo \(x)H(x - y, y) = 

x— >oo 

for each y and < l[ y oa )(x)H(x — y, y) < 1 for all x, y, we obtain 

/>oo 

lim f(x)u*(x) = lim / lr oo) (x)H(x - y, y)ip(y)u* {y)dy = 0, 
by the Lebesgue's dominated convergence theorem. Similarly, we conclude that 



lim 7(x)m*(x) = 0. 

X— 7-0 



We have 



o Jo 



r x rx—y 

H{x-y,y)p(y)u*(y)dy= I (p(y)u*(y)dy - \ \ h(z,y)dz(p(y)u*{y)dy. 
Thus, ju* G AC and 

(4.19) ^-{~f{x)u*{x)) = <p(x)u*(x) - / 
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The functions ipu* and v* arc intcgrablc. Consequently, if u* G L 1 then u* G T>(C) 
and C(u*) = 0. Since 

v*(x) = cp(x)u*(x) - 4-Mx)u*(x)) = -e- Q(3;) -^(e Q(x) 7(.T)u*(a;)), 
dx ax 

we obtain 

,*{x)e Qi - x) dx = - / — (e Q W~((x)u*(x))dx = e Q(z) -f(z)u* (z), 
Jz dx 

which shows that Kv*{y) = v*(y), by (|4.15|) . 

We now turn to the converse part. Suppose that u* is as in (|4.18[) . where v* is a 
fixed point of K . Since Q is non- increasing and v* is integrable, we see that 

lim j(x)u*(x) = 

x— >oo 

and that ipu* G L 1 . It is easily seen that u* satisfies equation (|4.19p . Integrating 
equation (|4. 1 9[) with respect to x from z to oo leads to 

OO J /"OO /"OO />X 

— (>y(x)u*(x))dx = / <p(x)u* (x)dx - / / h(x - y,y)(p(y)u*(y)dydx 
ax J z J z Jq 

and changing the order of integration in the last integral gives 

: pz poo 

h(x-y,y)tp(y)u*(y)dydx= / h(x - y,y)<p(y)u*(y)dxdy 



oo poo 



+ / h(x -y,y)ip(y)u*(y)dxdy. 

J z J y 

We have 

pz poo pz 

/ / h(x-y,y)tp(y)u*(y)dxdy = H(z - y,y)cp(y)u*(y)dy 

JO Jz Jo 

and 

fOO />OG /"OO 

/ K x -V^y) L P{y) u *{y)dxdy = j Lp(y)u*(y)dy. 

z J y J z 

Combining these we conclude that u* satisfies (|4.16|) . □ 

The following theorem guarantees that {P(t)}t>o is stochastic and its strong 
convergence to a unique stationary density u* that is given explicitly. 

Theorem 7. Suppose that the operator K as in (|4.15[) has an invariant density 
v* > a.e. and let 



1 



3 Qty)-Q(x) v *( y )dydx < OO. 



/o 70) J x 

Then the semigroup {P(t)}t>o is stochastic and for each initial density v we have 



where 



lim ||P(t)«-«*||i =0, 

t^-OO 



u *(x) = / e Q(v)-Ql*) v *( y )dy 

c r ){x) 
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is the unique stationary density of {P(t)} t > . 

Proof. By Proposition [5J the semigroup {P(t)}t>o is stochastic. From Theorem^ 
it follows that u* 6 T){C) and C(u*) = 0. Thus, u* is an invariant density for the 
stochastic semigroup {P{t)}t>o and it is unique, by Proposition [4] Since v*(x) > 
for a.e. x > 0, we conclude that u*(x) > for all x > 0. From assumptions (|4.3I) 
and (|4.4|) it follows that there is a 5o such that tp(y) > for y <E (0,<5o)- This and 
(147L7|) imply that 



OO POO 



p(x,y) ( p{y)dydx > 0, where p(x,y) = l( 0;X )(y)h(x - y,y). 

Consequently, we can find t > such that the operator P(t) is partially integral [14] 
and the result follows from Theorem [1] □ 

We conclude this section with sufficient conditions for mean ergodicity of the 
transition operator K. 

Proposition 8. Let K be a transition operator K with a bounded kernel k. 
Suppose that there exist a nonnegative measurable function V: (0, oo) — > [0, oo) which 
is bounded on bonded subsets of (0, oo) and constants a,d > such that 

/•oo 

(4.20) / V(x)k(x,y)dx<V(y)-l + al {0 , d) (y), y > 0. 

Then the operator K is mean ergodic on L 1 . 

Proof. Let Z n , n > 0, be a Markov chain with stochastic kernel /C given by 



K(y,B)= / k(x,y)dx, y > 0,B e £>((0, oo)). 
Jb 

Recall that a probability measure fi is invariant for the chain if and only if the measure 
/j, satisfies the equation 



H(B)= / K{y,B)n{dy) 
Jo 



for all Borel measurable sets B. We have 

fJ>i B ) = / / k(x,y)n(dy)dx. 
Jb Jo 

Thus each invariant probability measure is absolutely continuous with respect to the 
Lebesgue measure on (0, oo). Since K is the transition operator corresponding to /C, 
we have 



/•oo 

/ &{y,B)v(y)dy, BeB((0,oo)), 
Jo 



K 3 v{x)dx 
where K>{y,B) = fC(y, B) and 

K?{y,B)=\ ^-\z,B)1C{ y ,dz), y>0J>2. 



From Theorem 1 and Lemma 1 of |31| it follows that there exist a finite number of in- 
variant probability measures fx±, . . . , [An and a finite number of nonnegative functions 
Li, . . . Ljy such that 2i=i ^iiv) = 1 an( i 

n N 

(4.21) -J2K?(y,B)-^Y, L *(y)^ B ) 
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for all y and all Borel sets B. Let v\, . . . , «at be the densities of the invariant measures 
Hi,... ,Hn- Now let v S L 1 . From (|4.21[) and the Lebesgue dominated convergence 
theorem it follows that 

lim / - K j v(x)dx = / Li(y)v(y)dyvi(x)dx, 
n^oo J B n^-[ Jb^Jo 

for all Borel B. Moreover, the sequence ^ Y^j=i K-*v is bounded in L 1 . Thus, it is 
weakly convergent in L 1 and, by the mean ergodic theorem, it converges in L . □ 
We now apply the last result to our transition operator K. 

Corollary 9. Let K be the transition operator as in (|4.15|) with bounded h. 
Suppose that the function 



mi(y)= / xh(x,y)dx, y > 0, 
Jo 

is bounded on bounded subsets of (0, oo). If 

(4.22) limsupe Q ^ f ( mi (z)^4 - 1 I e~ Q(z) dz < 0, 

V^oo Jo V 7(2) / 

then the operator K is mean ergodic. 

Proof. Since K has kernel k given by (|4.9[) . we obtain 

k(x, y) < c x e Q{ ^ f ^rle- Q{z) dz =c u x, y > 0, 
Jo 7(2) 

where c\ is the upper bound for h. By Proposition [8] it is sufficient to check that 
the function V(x) = x, up to a multiplicative constant, satisfies condition (|4.20j) . We 
have 

/ V(x)h(x — z, z)dx = mi(z) + z, z > 0. 

Thus 

V(x)*(s, y)dx = e Qiy) [ (mi(«) + z)44 e_Q(2) ^ 
Jo 7(2) 



for all y > 0. Since 

2/ = 

we obtain 



0(v) [ V z VQ e -Q(*) dz + e Q(y) [ V e -Q(*)dz, 
Jo 7(2) 7o 



f °° V(x)fc(s, y)dr - V(y) = e Qfe) ( mUz)^- - l) e~ Q{z) dz, 
Jo V 7(2) / 



which is a bounded function on sets of the form (0,d). □ 
Remark 4- Observe that if Q(oo) = and 

mi(z)(p(z) 
lim sup — < 1 
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then condition (|4.22p holds, since we can find z a > and S > such that 

rmfz)^- -l<-5 for z > z , 
j(z) 



which implies that 

e Q(y) 



f TOl(z )<44 - l) e~ Q[z) dz < -ogOW-Odft.)^ _ yo ) 



for all y > yo > zq with the right-hand side going to — oo. 
If Q(oo) = — oo and 

limsup (mi (z) - ] < 

then condition (|4.22[) holds as well by d'Hospital's rule. 

4.2. Exponentially distributed bursts. Experimental findings in populations 
of cells indicate that the burst size is often exponentially distributed [33j so we now 
consider 

(4.23) h(x,y) = \e- x '\ x,y > 0, 

o 

where b > 0. The operator K as defined in (|4.15[) then takes the form 

Kv(x)= [ X l e ~(*~z)/b^l e -Q(z) [°° v *( y ) e Q(y) d ydz. 
Jo o 7(z) J z 

Note that the integrable function 

v *(x) = e -x/b-Q(x) 

is a fixed point of the operator K, since 

Kv*(x) = e-*/ h r 44e~ Q(2) ^ = e^/ 6 "^ . 
Jo 7(2) 

Again, an explicit stationary solution was recently obtained in [T], and we establish 
convergence to this stationary state with the following result. 

Corollary 10. Assume that conditions (|4.3[) and (|4.4[) hold and that h is 
exponential as in (|4.23p with b > 0. Suppose that 



(4.24) c := / — -e- :l;/b - Q(x) dx < 00, / eT x lh ~ Q{ ~ x) dx < 

Jo Jo 



00. 



Then the semigroup {P(t)} t >o is stochastic and for each initial density v we have 



lim \\P(t)v-u*\\i = 0, 

t—too 



where 



(4.25) u*(x) - -J— e -*/b-Q(x) 

cy{x) 
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is the unique stationary density of {P(t)} t > . 
Remark 5. Note that if Q(0) = oo and 

then the function x h-» e -' x / b -Q( x ) [ s integrablc on (0,oo). If, additionally, 

e -Q(x) rS 

liminf7(x) > 0, lim — — — < oo, and / j(x) r dx < oo 
x-^oo x^o ^(xy J Q 

for some S, r > 0, then condition (|4.24[) holds. Furthermore, if it should happen that 
b, j(x) and u*{x) are known or can be approximated from data, then it is possible to 
estimate tp(x) from 

. . 1 , . (^f(x)u* (x))' 
(4.26) if{x) = - 7 (.t) + m ^ ;; ■ 

Finally, note that if tp is assumed to be bounded, then u* has an exponential tail, 
from which we can deduce the parameter b. 

Remark 6 (Bifurcation in the continuous case/ As in the discrete formulation of 
the model we can use relation (|4.26|) to derive bifurcation properties of the stationary 
density as a function of the relevant parameters. Namely, the number of extrema are 
linked to the number of solutions of 

( \ ! U \ 

<PW = —r- +7 W- 



In all examples below, we consider a linear degradation function, 7(2;) = 72; with 
7 > 0. 

Example 6. Consider the function (p of the form 

/ s x l + ©s* A / A\ 1 
(4-27) <P(x) = a , ^ N = X + A 1 - x 



A + Ax N A V A/ A + Ax N ' 
where A, A, A, N are positive constants and > 0. Then 

Q(x) = c - A l 0g{x) + bg(A + Ax N ), 

where c\ is a constant. The stationary density is given by 
(4.28) u*{x) = (c-/)- 1 e- x / b x x ^ A ^ 1 - 1 (A + Ax N ) e , 

where 

( 4 - 29 ) = 777- (© - V 

v ; AA7 V A 

This solution has been extensively studied in terms of numbers of maxima (P-bifurcation)| 
in [15] when = 1. When = A = A = 1 the density u* is that of a gamma distri- 
bution, as obtained in [5J. 

Example 7. Consider the case of linear regulation with the function ip of the form 

<p(x) — A + \x, 
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where Ao, A arc nonnegativc constants. If 

(4.30) T > - and A > 0, 

o 7 

then u* is intcgrable and is given by the gamma distribution 

<«» «*<*> = ^(HH^^ 

which is a continuous approximation of the negative binomial distribution previously 
obtained, as in [27] . 

4.3. Other examples. In this subsection we consider some more exactly solv- 
able examples. The class of examples we provide generalizes the exponentially dis- 
tributed case of h. Let v{y) be a positive, decreasing, and absolutely continuous 
function on (0, oo) such that u(y) — > as y — > oo. Consider the function 

(4.32) h{x,y) = - V ' {X + V \ y,z>0. 

Then for each y the function x t— > h(x, y) is a density and 

ny) 

The operator K can be thus rewritten as 

Kv(x) = - f X ^^j e-QW [°° v(y)e^dydz. 

JO v \ z ) l\ z ) Jz 

It is easily seen that if the function 

v *(x) = - V '(x)e- Q{x) 

is integrable then Kv*{x) = v*(x), thus we obtain the following. 
Corollary 11. Let h be as in (|4.32|) . Suppose that 

f°° v(x) f°° 
c := — Y-{e~ Q ^dx < oo and - / i/'(x)e~ Q ^dx < oo. 
Jo 7W 7o 

T/ien i/ie semigroup {P(i)} t > is stochastic and for each initial density v we have 



where 



lim ||P(t)u-«*||i =0, 



= — ---e 

cy(x) 

is the unique stationary density of {P(t)} t >o- 

Remark 7 (^Bifurcation in the continuous case-again,). As before (see Remark |6]), 
the number of extrema are linked to the number of solutions of 

v{x) 
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Note that if it should happen that v(x),j(x) and u*(x) are known or can be approx- 
imated from data, then it is possible to estimate p{x) from 

uyx) u*(x) 

If rrii(x) — zh(z,x)dz < oo, then only the knowledge of m\{x) is sufficient as 

v'(x) _ 1 + m[(x) 
u(x) mi(x) 

In the examples below, we take a linear degradation function j(x) = 72;, with 7 > 0. 
Example 8. Suppose that the function v is of the form 

v(x) = (a + x)~^ 

where a, (3 > and that the function ip is of the form (|4.27j) . If 

n ^© 

P> -T- + 1 
7A 

then the assumptions of Corollary [TT] are satisfied and the stationary density u* is 
given by 

where 9 is as in (|4.29J) . 

Example 9. Suppose that the function v is of the form 

v (x) = e ~(^+^ 2 ) ; 

where a, f3 > 0. Consider the case of linear regulation with the function tp of the form 

p(x) = A + Xix, 
where Ao, Ai are nonnegative constants. If 

A > 0, 

then u* is integrable and is given by 

u (X) = — x t e v y 
07 

Consider the case of quadratic regulation with the function ip of the form 

<p(x) = A + Xix + X 2 x 2 , 
where Ao, Ai, A2 are nonnegative constants. If 

A 2 

P > — and A > 0, 
2 7 
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then u* is integrable and is given by 

u * {x) = J_ x ^-i e -(-^)-(^)- 2 . 

07 



Example 10. Suppose that the function v is of the form 

v(x) = (a — x)P , 

where a,/3 > 0, for all x < a, and fix) = for x > a. Suppose the function ip is given 
by (|4.27p where A, A, A, AT are positive constants and > 0. Then the stationary 
density u* is integrable and is given by, for all x < a, 

u*(x) = —(a - xfx^V 1 - 1 ^ + Ax N ) e , 
C7 

where 9 is as in (|4.29|) . Convergence is obtained in the state space (0,a). 

5. Conclusions and summary. In this paper we have presented both a discrete 
Markov process formulation as well as a continuous model formulation for bursting 
gene expression. Our development of the discrete model formulation in Section 13.11 
allowed us to prove a very general convergence result in Theorem [2] and then to use 
that result to explore a variety of examples in Section 13.21 when the burst amplitude 
is geometrically distributed. In Section @] we developed the analogous continuous 
model for bursting expression. Section 14.11 contains the general development with 
Proposition |4] limiting the number of invariant densities of the semigroup {P(t)} t >o, 
while Proposition [5] uses mean ergodicity of the transition operator K to show that 
{P(t)}t>o is stochastic. Theorems |6] and [7] give criteria for a unique stationary den- 
sity u* of {P(i)}t>o and for convergence to that stationary density. In Section |4~21 
we have used these results in a number of specific examples when the burst ampli- 
tudes are exponentially distributed-a situation often noted experimentally. Section 
14.31 concludes with an examination of a generalization of the exponential distribution 
of burst amplitudes. 
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